## REPLICATION 
## Business Owners and Executives as Politicians: The Effect on Public Policy
## produces & saves data frames of RD results 
## includes validity tests and robustness checks (located in Supplemental Information)


rm(list=ls(all=TRUE))


library(rdrobust)
library(stargazer)
library(tidyverse)
library(foreign)
library(rdd)
library(xtable)

## set working directory
setwd("")

source("rddensity_fun.R")
source("rdbwdensity.R")
source("rddensity.R")
source("replication_mayors_source.R")


# mayoral & fiscal data 
load("replication_data_business_mayors_RD.RData")

## subset to elections with one business executive
mayors_subs <- subset(mayors, !is.na(exec_margin))


#  Density test 
rddensity(X = mayors_subs$exec_margin, vce="jackknife")

#### results
# Method                        T                   P > |T|             
# Robust Bias-Corrected         -1.3196             0.187 

#### McCrary Test
DCdensity(runvar = mayors_subs$exec_margin
          , cutpoint = 0
          , bin = NULL
          , bw = NULL
          , verbose = TRUE
          , plot = TRUE
          , ext.out = FALSE )

#### results 
# Using calculated bin size:  0.010 
# Using calculated bandwidth:  0.119 
# Log difference in heights is  -0.187  with SE  0.202 
# this gives a z-stat of  -0.925 
# and a p value of  0.355 
# [1] 0.3550962

#########################################
#####  CONTINUITY & PLACEBO TESTS  ######
#########################################


covlist <- c("pc_ttl_revs_lag1", "pc_ttl_exps_lag1", "pc_ttl_own_revs_lag1",
             "pc_ttl_taxes_lag1", "pc_ttl_sales_taxes_lag1", "pc_ttl_prop_taxes_lag1",
             "pc_ttl_chgs_miscrev_lag1",
             "pc_ttl_debt_lag1", "pc_ttl_debt_issued_lag1", "pc_st_debt_lag1",
             "pc_police_lag1", "pc_fire_lag1", "pc_admin_lag1", "pc_waste_lag1",
             "pc_roads_lag1", "pc_parks_lag1", "pc_libraries_lag1",
             "pc_health_lag1", "pc_housing_lag1","pc_welfare_lag1",
             "form_alt",
             "ipctwhite", "iunemp", "ihomeown", "imedhhinc")




results_balance <- data.frame()

### covariate-adjusted placebo tests


covlist_fiscal <- c("pc_ttl_revs", "pc_ttl_exps", "pc_ttl_own_revs",
                    "pc_ttl_taxes", "pc_ttl_sales_taxes", "pc_ttl_prop_taxes",
                    "pc_ttl_chgs_miscrev", 
                    "pc_ttl_debt", "pc_ttl_debt_issued", "pc_st_debt",
                    "pc_police", "pc_fire", "pc_admin", "pc_waste", 
                    "pc_roads", "pc_parks", "pc_libraries",
                    "pc_health", "pc_housing","pc_welfare")



h5 <- 0.05  # 5% bandwidth
wgts <- (1 - abs(mayors_subs$exec_margin / h5)) * (abs(mayors_subs$exec_margin) <= h5)


models_covs <- lapply(covlist_fiscal, function(x) {
  

  lm(substitute(i ~ (exec_margin > 0) * exec_margin + 
                  j +
                  ipopulation_log_lag2 +
                  ipctwhite_lag2 +
                  imedhhinc_log_lag2 +
                  imedhouseval_log_lag2,
                list(i = as.name(paste0(x, "_lag1")), 
                     j = as.name(paste0(x, "_lag4"))
                     )), 
     data = mayors_subs, weights = wgts, subset = wgts > 0)
  
})

covs_conv <- lapply(models_covs, summary)

covs_robust <- lapply(models_covs, function(x) {
  
  coeftest(x, vcov=vcovHC(x, type="HC1")) 
})



for(x in 1:length(covlist_fiscal)){
  temp_results <- data.frame(cbind(as.character(paste0(covlist_fiscal[[x]], "_lag1")), 
                                   format_num(covs_conv[[x]]$coefficients[2,1]),
                                   format_num(covs_robust[[x]][2,2]),
                                   format_num(covs_robust[[x]][2,4]),
                                   format_num(h5, digits = 3),
                                   format_num(covs_conv[[x]]$df[2] + covs_conv[[x]]$df[3], digits = 0),
                                   "covs"))

  names(temp_results) <- c("outcome", "coef", "robust_se", "robust_pval", "h", "n", "specs")
  
  results_balance <- rbind(results_balance, temp_results)
}








rdrobust_bal_specs_covs <- lapply(covlist_fiscal, function(x) {
  
  substitute(with(mayors_subs, rdrobust(i, exec_margin,
                      covs = cbind(j, 
                                   ipopulation_log_lag2,
                                   ipctwhite_lag2,
                                   imedhhinc_log_lag2,
                                   imedhouseval_log_lag1))),
                                   list(i = as.name(paste0(x, "_lag1")), 
                                        j = as.name(paste0(x, "_lag4"))
                                   ))
})



rdrobust_bal_tests_covs <- lapply(rdrobust_bal_specs_covs, eval)


## these results use CCT optimal bandwidth but regular SEs
for(j in 1:length(covlist_fiscal)){
  
  ## covs model    
  h_CCT <- as.numeric(rdrobust_bal_tests_covs[[j]]$bws[1, 2])
  
  mayors_subs$wgts_CCT <- weights_CCT(mayors_subs$exec_margin, h_CCT)
  
  
  covs <-    lm(substitute(i ~ (exec_margin > 0) * exec_margin +
                             j +
                             ipopulation_log_lag2 +
                             ipctwhite_lag2 +
                             imedhhinc_log_lag2 +
                             imedhouseval_log_lag1,
                           list(i = as.name(paste0(covlist_fiscal[j], "_lag1")), 
                                j = as.name(paste0(covlist_fiscal[j], "_lag4"))
                           )), 
                     data = mayors_subs, weights = wgts_CCT, subset = wgts_CCT > 0)
  
  summ_covs <- summary(covs)
  
  covs_robust <- coeftest(covs, vcov=vcovHC(covs, type="HC1")) 
  

  temp_results <- data.frame(cbind(as.character(covlist[j]), 
                                   format_num(summ_covs$coefficients[2,1]),
                                   format_num(covs_robust[2,2]),
                                   format_num(covs_robust[2,4]),
                                   format_num(h_CCT),
                                   format_num(summ_covs$df[2] + summ_covs$df[3], digits = 0),
                                   "covs"))

  names(temp_results) <- c("outcome", "coef", "robust_se", "robust_pval", "h", "n", "specs")
  
  results_balance <- rbind(results_balance, temp_results)
  
}




## placebo tests for city characteristics w/ covariates - 5%

covlist_ses <-  c("population", "pctwhite", "unemp", "homeown", "medhhinc", "medhouseval", "form_alt")

h5 <- 0.05  # 5% bandwidth
mayors_subs$wgts5 <- (1 - abs(mayors_subs$exec_margin / h5)) * (abs(mayors_subs$exec_margin) <= h5)

population_lm <- lm(ipopulation_log ~ (exec_margin > 0) * exec_margin +
                      ipopulation_log_lag3 +
                      ipctwhite +
                      imedhhinc_log +
                      imedhouseval_log,
                    data = mayors_subs, weights = wgts5, subset = wgts5 > 0
)

                   

pctwhite_lm <- lm(ipctwhite ~ (exec_margin > 0) * exec_margin +
                      ipopulation_log +
                      ipctwhite_lag3 +
                      imedhhinc_log +
                      imedhouseval_log,
                    data = mayors_subs, weights = wgts5, subset = wgts5 > 0
)


medhhinc_lm <- lm(imedhhinc_log ~ (exec_margin > 0) * exec_margin +
                    ipopulation_log +
                    ipctwhite +
                    imedhhinc_log_lag3 +
                    imedhouseval_log,
                  data = mayors_subs, weights = wgts5, subset = wgts5 > 0
)



medhouseval_lm <- lm(imedhouseval_log ~ (exec_margin > 0) * exec_margin +
                        ipopulation_log +
                        ipctwhite +
                        imedhhinc_log +
                        imedhouseval_log_lag3,
                      data = mayors_subs, weights = wgts5, subset = wgts5 > 0
)



unemp_lm <- lm(iunemp ~ (exec_margin > 0) * exec_margin +
                        ipopulation_log +
                        ipctwhite +
                        imedhhinc_log +
                        imedhouseval_log +
                        iunemp_lag3,
                      data = mayors_subs, weights = wgts5, subset = wgts5 > 0
)


homeown_lm <- lm(ihomeown ~ (exec_margin > 0) * exec_margin +
                  ipopulation_log +
                  ipctwhite +
                  imedhhinc_log +
                  imedhouseval_log +
                  ihomeown_lag3,
                data = mayors_subs, weights = wgts5, subset = wgts5 > 0
)



form_alt_lm <- lm(form_alt ~ (exec_margin > 0) * exec_margin +
                 ipopulation_log +
                 ipctwhite +
                 imedhhinc_log +
                 imedhouseval_log +
                 form_alt_lag3,
               data = mayors_subs, weights = wgts5, subset = wgts5 > 0
)




for(i in 1:length(covlist_ses)){
  model <- eval(as.name(paste0(covlist_ses[i], "_lm")))
  
  temp_results <- data.frame(cbind(
    covlist_ses[i],
    format_num(summary(model)$coefficients[2, 1], digits = 4), # coef
    format_num(coeftest(model, vcov=vcovHC(model, type="HC1"))[2, 2], digits = 3), # se
    format_num(coeftest(model, vcov=vcovHC(model, type="HC1"))[2, 4], digits = 3), # pval
    format_num(h5, digits = 3), # h
    summary(model)$df[1] + summary(model)$df[2], # n
    "covs"))

  
  names(temp_results) <- c("outcome", "coef", "robust_se", "robust_pval", "h", "n", "specs")
  
  results_balance <- rbind(results_balance, temp_results)
  
}




####  covariates with CCT optimal bandwidth

population_rdrobust <- with(mayors_subs, rdrobust(ipopulation_log, exec_margin,
                                                  covs = cbind(ipopulation_log_lag3,
                                                               ipctwhite,
                                                               imedhhinc_log,
                                                               imedhouseval_log)
                                                  ))

h_CCT_population <- population_rdrobust$bws[1, 2]

mayors_subs$wgts_CCT_population <- (1 - abs(mayors_subs$exec_margin / h_CCT_population)) * (abs(mayors_subs$exec_margin) <= h_CCT_population)

population_lm <- lm(ipopulation_log ~ (exec_margin > 0) * exec_margin +
                      ipopulation_log_lag3 +
                      ipctwhite +
                      imedhhinc_log +
                      imedhouseval_log,
                    data = mayors_subs, weights = wgts_CCT_population, subset = wgts_CCT_population > 0
)


pctwhite_rdrobust <- with(mayors_subs, rdrobust(ipctwhite, exec_margin,
                                             covs = cbind(ipopulation_log,
                                                          ipctwhite_lag3,
                                                          imedhhinc_log,
                                                          imedhouseval_log)
                                             ))

h_CCT_pctwhite <- pctwhite_rdrobust$bws[1, 2]

mayors_subs$wgts_CCT_pctwhite <- (1 - abs(mayors_subs$exec_margin / h_CCT_pctwhite)) * (abs(mayors_subs$exec_margin) <= h_CCT_pctwhite)

pctwhite_lm <- lm(ipctwhite ~ (exec_margin > 0) * exec_margin +
                    ipopulation_log +
                    ipctwhite_lag3 +
                    imedhhinc_log +
                    imedhouseval_log,
                  data = mayors_subs, weights = wgts_CCT_pctwhite, subset = wgts_CCT_pctwhite > 0
)





medhhinc_rdrobust <- with(mayors_subs, rdrobust(imedhhinc_log, exec_margin,
                                             covs = cbind(ipopulation_log,
                                                          ipctwhite,
                                                          imedhhinc_log_lag3,
                                                          imedhouseval_log)
                                             ))

h_CCT_medhhinc <- medhhinc_rdrobust$bws[1, 2]

mayors_subs$wgts_CCT_medhhinc <- (1 - abs(mayors_subs$exec_margin / h_CCT_medhhinc)) * (abs(mayors_subs$exec_margin) <= h_CCT_medhhinc)



medhhinc_lm <- lm(imedhhinc_log ~ (exec_margin > 0) * exec_margin +
                    ipopulation_log +
                    ipctwhite +
                    imedhhinc_log_lag3 +
                    imedhouseval_log,
                  data = mayors_subs, weights = wgts_CCT_medhhinc, subset = wgts_CCT_medhhinc > 0
)




medhouseval_rdrobust <- with(mayors_subs, rdrobust(imedhouseval_log, exec_margin,
                                             covs = cbind(ipopulation_log,
                                                          ipctwhite,
                                                          imedhhinc_log,
                                                          imedhouseval_log_lag3)
                                             ))

h_CCT_medhouseval <- medhouseval_rdrobust$bws[1, 2]

mayors_subs$wgts_CCT_medhouseval <- (1 - abs(mayors_subs$exec_margin / h_CCT_medhouseval)) * (abs(mayors_subs$exec_margin) <= h_CCT_medhouseval)


medhouseval_lm <- lm(imedhouseval_log ~ (exec_margin > 0) * exec_margin +
                       ipopulation_log +
                       ipctwhite +
                       imedhhinc_log +
                       imedhouseval_log_lag3,
                     data = mayors_subs, weights = wgts_CCT_medhouseval, subset = wgts_CCT_medhouseval > 0
)



unemp_rdrobust <- with(mayors_subs, rdrobust(iunemp, exec_margin,
                                                covs = cbind(ipopulation_log,
                                                             ipctwhite,
                                                             imedhhinc_log,
                                                             imedhouseval_log,
                                                             iunemp_lag3)
                                          ))
h_CCT_unemp <- unemp_rdrobust$bws[1, 2]


mayors_subs$wgts_CCT_unemp <- (1 - abs(mayors_subs$exec_margin / h_CCT_unemp)) * (abs(mayors_subs$exec_margin) <= h_CCT_unemp)



unemp_lm <- lm(iunemp ~ (exec_margin > 0) * exec_margin +
                 ipopulation_log +
                 ipctwhite +
                 imedhhinc_log +
                 imedhouseval_log +
                 iunemp_lag3,
               data = mayors_subs, weights = wgts_CCT_unemp, subset = wgts_CCT_unemp > 0
)




homeown_rdrobust <- with(mayors_subs, rdrobust(ihomeown, exec_margin,
                                                covs = cbind(ipopulation_log,
                                                             ipctwhite,
                                                             imedhhinc_log,
                                                             imedhouseval_log,
                                                             ihomeown_lag3)
                                            ))

h_CCT_homeown <- homeown_rdrobust$bws[1, 2]

mayors_subs$wgts_CCT_homeown <- (1 - abs(mayors_subs$exec_margin / h_CCT_homeown)) * (abs(mayors_subs$exec_margin) <= h_CCT_homeown)


homeown_lm <- lm(ihomeown ~ (exec_margin > 0) * exec_margin +
                   ipopulation_log +
                   ipctwhite +
                   imedhhinc_log +
                   imedhouseval_log +
                   ihomeown_lag3,
                 data = mayors_subs, weights = wgts_CCT_homeown, subset = wgts_CCT_homeown > 0
)




form_alt_rdrobust <- with(mayors_subs, rdrobust(form_alt, exec_margin,
                                            covs = cbind(ipopulation_log,
                                                         ipctwhite,
                                                         imedhhinc_log,
                                                         imedhouseval_log,
                                                         form_alt_lag3)
                                         ))

h_CCT_form_alt <- form_alt_rdrobust$bws[1, 2]

mayors_subs$wgts_CCT_form_alt <- (1 - abs(mayors_subs$exec_margin / h_CCT_form_alt)) * (abs(mayors_subs$exec_margin) <= h_CCT_form_alt)

form_alt_lm <- lm(form_alt ~ (exec_margin > 0) * exec_margin +
                    ipopulation_log +
                    ipctwhite +
                    imedhhinc_log +
                    imedhouseval_log +
                    form_alt_lag3,
                  data = mayors_subs, weights = wgts_CCT_form_alt, subset = wgts_CCT_form_alt > 0
)





for(i in 1:length(covlist_ses)){
  model <- eval(as.name(paste0(covlist_ses[i], "_lm")))
  wgts_CCT_var <- as.name(paste0("mayors_subs$wgts_CCT_", covlist_ses[i]))
  mayors_subs$wgts_CCT <- eval(parse(text = wgts_CCT_var))
  
  temp_results <- data.frame(cbind(
    covlist_ses[i],
    format_num(summary(model)$coefficients[2, 1], digits = 4), # coef
    format_num(coeftest(model, vcov=vcovHC(model, type="HC1"))[2, 2], digits = 3), # se
    format_num(coeftest(model, vcov=vcovHC(model, type="HC1"))[2, 4], digits = 3), # pval
    format_num(eval(as.name(paste0("h_CCT_", covlist_ses[i]))), digits = 3), # h
    summary(model)$df[1] + summary(model)$df[2], # n
    "covs"))
  
  names(temp_results) <- c("outcome", "coef", "robust_se", "robust_pval", "h", "n", "specs")
  
  results_balance <- rbind(results_balance, temp_results)
  
}


## results - optimal bandwidth
results_balance_covs_CCT <- results_balance[c(21:40, 48:54),]
## results - 5% bandwidth
results_balance_covs_05 <- results_balance[c(1:20, 41:47),]

save(list = c("results_balance_covs_CCT", "results_balance_covs_05"), file = "replication_results_covs_balance_tests.RData")



## continuity tests using rdrobust
results_balance_rdrobust <- data.frame()


## fiscal covs
for (x in 1:length(covlist_fiscal)) {
  
  temp_results <- data.frame(cbind(as.character(covlist_fiscal[[x]]), 
                                   format_num(rdrobust_bal_tests_covs[[x]]$coef[1]),
                                   format_num(rdrobust_bal_tests_covs[[x]]$pv[3]),
                                   format_num(rdrobust_bal_tests_covs[[x]]$ci[3,1]),
                                   format_num(rdrobust_bal_tests_covs[[x]]$ci[3,2]),
                                   format_num(rdrobust_bal_tests_covs[[x]]$bws[1, 2]),
                                   format_num(rdrobust_bal_tests_covs[[x]]$Nh[1] + rdrobust_bal_tests_covs[[x]]$Nh[2], digits = 0),
                                   "covs"))
  
  names(temp_results) <- c("outcome", "coef", "p_val", "ci_low", "ci_upp", 
                           "h", "n", "specs")
  
  results_balance_rdrobust <- rbind(results_balance_rdrobust, temp_results)
  
}





for(i in 1:length(covlist_ses)){
  rdrobust_model <- eval(as.name(paste0(covlist_ses[i], "_rdrobust")))



  temp_results <- data.frame(cbind(as.character(covlist_ses[i]),
                                   format_num(rdrobust_model$coef[1]),
                                   format_num(rdrobust_model$pv[3]),
                                   format_num(rdrobust_model$ci[3,1]),
                                   format_num(rdrobust_model$ci[3,2]),
                                   format_num(rdrobust_model$bws[1, 2]),
                                   format_num(rdrobust_model$Nh[1] + rdrobust_model$Nh[2], digits = 0),
                                   "covs"))

  names(temp_results) <- c("outcome", "coef", "p_val", "ci_low", "ci_upp",
                           "h", "n", "specs")

  results_balance_rdrobust <- rbind(results_balance_rdrobust, temp_results)

}


results_balance_rdrobust_covs <- filter(results_balance_rdrobust, specs == "covs")

save(results_balance_rdrobust_covs, file = "replication_results_covs_balance_tests_rdrobust.RData")



##########################################
###### QUADRATIC RD SPECIFICATIONS  ######
##########################################

## List of DVs
varlist <- c("ttl_revs", "ttl_own_revs", 
             "ttl_taxes", "ttl_sales_taxes", "ttl_prop_taxes",
             "ttl_chgs_miscrev", 
             "ttl_debt", "ttl_debt_issued", "st_debt",
             "ttl_exps", 
             "police", "fire", "admin", "waste", 
             "roads", "parks", "libraries",
             "health", "housing","welfare",
             "basic_services", "developmental", "redistributive")





############ CCT bandwidths w/ quadratics 

results_CCT <- data.frame()


## w/ covs

rdrobust_specs_covs <- lapply(varlist, function(x) {
  
  substitute(with(mayors_subs, rdrobust(i, exec_margin,
                                        covs = cbind(j, 
                                                     log(ipopulation),
                                                     ipctwhite,
                                                     log(imedhhinc),
                                                     log(imedhouseval)), p = 2)),
             list(i = as.name(paste0("pc_", x, "_lead2")), 
                  j = as.name(paste0("pc_", x, "_lag1"))
             ))
})



rdrobust_results_covs <- lapply(rdrobust_specs_covs, eval)



## these results use CCT optimal bandwidth but regular SEs
for(j in 1:length(varlist)){
  
  ## covs model    
  h_CCT <- as.numeric(rdrobust_results_covs[[j]]$bws[1, 2])
  
  mayors_subs$wgts_CCT <- weights_CCT(mayors_subs$exec_margin, h_CCT)
  
  
  covs <-    lm(substitute(i ~ (exec_margin > 0) * exec_margin + (exec_margin > 0) * I(exec_margin^2) +
                             j +
                             log(ipopulation) +
                             ipctwhite +
                             log(imedhhinc) +
                             log(imedhouseval),
                           list(i = as.name(paste0("pc_", varlist[j], "_lead2")), 
                                j = as.name(paste0("pc_", varlist[j], "_lag1"))
                           )), 
                data = mayors_subs, weights = wgts_CCT, subset = wgts_CCT > 0)
  
  summ_covs <- summary(covs)
  
  covs_robust <- coeftest(covs, vcov=vcovHC(covs, type="HC1")) 
  
  
  temp_results <- data.frame(cbind(as.character(varlist[j]), 
                                   format_num(summ_covs$coefficients[2,1]),
                                   format_num(covs_robust[2,2]),
                                   format_num(covs_robust[2,4]),
                                   format_num(h_CCT),
                                   format_num(summ_covs$df[2] + summ_covs$df[3], digits = 0),
                                   "CCT covs"))

  names(temp_results) <- c("outcome", "coef", "robust_se", "robust_pval", "h", "n", "specs")
  
  results_CCT <- rbind(results_CCT, temp_results)
  
}




results_CCT_quadratic <- results_CCT


save(results_CCT_quadratic, file = "replication_results_CCT_quadratic_specs.RData")





##################################################
###### ALTERNATIVE RD SPECIFICATIONS W/     ######
###### BIAS-CORRECTED CONFIDENCE INTERVALS  ######
##################################################

#### Estimates with bias-corrected confidence intervals (using RDrobust package)


## List of DVs
varlist <- c("ttl_revs", "ttl_own_revs", 
             "ttl_taxes", "ttl_sales_taxes", "ttl_prop_taxes",
             "ttl_chgs_miscrev", 
             "ttl_debt", "ttl_debt_issued", "st_debt",
             "ttl_exps", 
             "police", "fire", "admin", "waste", 
             "roads", "parks", "libraries",
             "health", "housing","welfare",
             "basic_services", "developmental", "redistributive")

dv_means <- lapply(varlist, function(x) {

  substitute(mean(mayors_subs$i, na.rm=TRUE), list(i = as.name(paste0("pc_", x, "_lead2"))))

})

dv_std_dev <- lapply(varlist, function(x) {

  substitute(sd(mayors_subs$i, na.rm=TRUE), list(i = as.name(paste0("pc_", x, "_lead2"))))

})


models_base <- lapply(varlist, function(x) {

  substitute(rdrobust(mayors_subs$i, mayors_subs$exec_margin), list(i = as.name(paste0("pc_", x, "_lead2"))))

})


models_log <- lapply(varlist, function(x) {

  substitute(rdrobust(log(mayors_subs$i + 1), mayors_subs$exec_margin), list(i = as.name(paste0("pc_", x, "_lead2"))))

})


models_covs <- lapply(varlist, function(x) {
  
  substitute(rdrobust(mayors_subs$i, mayors_subs$exec_margin, 
                      covs = cbind(mayors_subs$j, log(mayors_subs$ipopulation),
                                   mayors_subs$ipctwhite, log(mayors_subs$imedhhinc), 
                                   log(mayors_subs$imedhouseval))),
             list(i = as.name(paste0("pc_", x, "_lead2")), j = as.name(paste0("pc_", x, "_lag1"))))       
  
  
})



models_log_covs <- lapply(varlist, function(x) {

  substitute(rdrobust(log(mayors_subs$i + 1), mayors_subs$exec_margin,
                      covs = cbind(log(mayors_subs$j +1), log(mayors_subs$ipopulation),
                                   mayors_subs$ipctwhite, log(mayors_subs$imedhhinc),
                                   log(mayors_subs$imedhouseval))),
             list(i = as.name(paste0("pc_", x, "_lead2")), j = as.name(paste0("pc_", x, "_lag1"))))


})


models_diff <- lapply(varlist, function(x) {

  substitute(rdrobust((mayors_subs$i - mayors_subs$j), mayors_subs$exec_margin),
             list(i = as.name(paste0("pc_", x, "_lead2")), j = as.name(paste0("pc_", x, "_lag1"))))


})



models_diff_covs <- lapply(varlist, function(x) {

  substitute(rdrobust((mayors_subs$i - mayors_subs$j),  mayors_subs$exec_margin,
                      covs = cbind(log(mayors_subs$ipopulation),
                                   mayors_subs$ipctwhite, log(mayors_subs$imedhhinc),
                                   log(mayors_subs$imedhouseval))),
             list(i = as.name(paste0("pc_", x, "_lead2")), j = as.name(paste0("pc_", x, "_lag1"))))


})


outcome_means <- lapply(dv_means, eval)
outcome_std_dev <- lapply(dv_std_dev, eval)

rd_ests_base <- lapply(models_base, eval)

rd_ests_log <- lapply(models_log, eval)

rd_ests_log_covs <- lapply(models_log_covs, eval)

rd_ests_covs <- lapply(models_covs, eval)

rd_ests_diff <- lapply(models_diff, eval)

rd_ests_diff_covs <- lapply(models_diff_covs, eval)






results <- data.frame()

for (x in 1:length(varlist)) {
  
  
  temp_results <- data.frame(cbind(as.character(varlist[[x]]),
                                   format_num(rd_ests_base[[x]]$coef[1]),
                                   format_num(rd_ests_base[[x]]$pv[3]),
                                   format_num(rd_ests_base[[x]]$ci[3,1]),
                                   format_num(rd_ests_base[[x]]$ci[3,2]),
                                   format_num(rd_ests_base[[x]]$bws[1, 2]),
                                   format_num(rd_ests_base[[x]]$Nh[1] + rd_ests_base[[x]]$Nh[2], digits = 0),
                                   format_num(outcome_means[[x]]),
                                   format_num(outcome_std_dev[[x]]),
                                   "base"))
  names(temp_results) <- c("outcome", "coef", "p_val", "ci_low", "ci_upp",
                           "h", "n", "mean", "sd", "specs")

  results <- rbind(results, temp_results)



  temp_results <- data.frame(cbind(as.character(varlist[[x]]),
                                   format_num(rd_ests_log[[x]]$coef[1]),
                                   format_num(rd_ests_log[[x]]$pv[3]),
                                   format_num(rd_ests_log[[x]]$ci[3,1]),
                                   format_num(rd_ests_log[[x]]$ci[3,2]),
                                   format_num(rd_ests_log[[x]]$bws[1, 2]),
                                   format_num(rd_ests_log[[x]]$Nh[1] + rd_ests_log[[x]]$Nh[2], digits = 0),
                                   format_num(outcome_means[[x]]),
                                   format_num(outcome_std_dev[[x]]),
                                   "log"))

  names(temp_results) <- c("outcome", "coef", "p_val", "ci_low", "ci_upp",
                           "h", "n", "mean", "sd", "specs")

  results <- rbind(results, temp_results)

  temp_results <- data.frame(cbind(as.character(varlist[[x]]),  
                                   format_num(rd_ests_covs[[x]]$coef[1]),
                                   format_num(rd_ests_covs[[x]]$pv[3]),
                                   format_num(rd_ests_covs[[x]]$ci[3,1]),
                                   format_num(rd_ests_covs[[x]]$ci[3,2]),
                                   format_num(rd_ests_covs[[x]]$bws[1, 2]),
                                   format_num(rd_ests_covs[[x]]$Nh[1] + rd_ests_covs[[x]]$Nh[2], digits = 0),
                                   format_num(outcome_means[[x]]),
                                   format_num(outcome_std_dev[[x]]),
                                   "covs"))
  
  names(temp_results) <- c("outcome", "coef", "p_val", "ci_low", "ci_upp", 
                           "h", "n", "mean", "sd", "specs")
  
  results <- rbind(results, temp_results)
  
  temp_results <- data.frame(cbind(as.character(varlist[[x]]),
                                   format_num(rd_ests_log_covs[[x]]$coef[1]),
                                   format_num(rd_ests_log_covs[[x]]$pv[3]),
                                   format_num(rd_ests_log_covs[[x]]$ci[3,1]),
                                   format_num(rd_ests_log_covs[[x]]$ci[3,2]),
                                   format_num(rd_ests_log_covs[[x]]$bws[1, 2]),
                                   format_num(rd_ests_log_covs[[x]]$Nh[1] + rd_ests_log_covs[[x]]$Nh[2], digits = 0),
                                   format_num(outcome_means[[x]]),
                                   format_num(outcome_std_dev[[x]]),
                                   "log covs"))

  names(temp_results) <- c("outcome", "coef", "p_val", "ci_low", "ci_upp",
                           "h", "n", "mean", "sd", "specs")

  results <- rbind(results, temp_results)

  temp_results <- data.frame(cbind(as.character(varlist[[x]]),
                                   format_num(rd_ests_diff[[x]]$coef[1]),
                                   format_num(rd_ests_diff[[x]]$pv[3]),
                                   format_num(rd_ests_diff[[x]]$ci[3,1]),
                                   format_num(rd_ests_diff[[x]]$ci[3,2]),
                                   format_num(rd_ests_diff[[x]]$bws[1, 2]),
                                   format_num(rd_ests_diff[[x]]$Nh[1] + rd_ests_diff[[x]]$Nh[2], digits = 0),
                                   format_num(outcome_means[[x]]),
                                   format_num(outcome_std_dev[[x]]),
                                   "differenced dv"))

  names(temp_results) <- c("outcome", "coef", "p_val", "ci_low", "ci_upp",
                           "h", "n", "mean", "sd", "specs")

  results <- rbind(results, temp_results)

  temp_results <- data.frame(cbind(as.character(varlist[[x]]),
                                   format_num(rd_ests_diff_covs[[x]]$coef[1]),
                                   format_num(rd_ests_diff_covs[[x]]$pv[3]),
                                   format_num(rd_ests_diff_covs[[x]]$ci[3,1]),
                                   format_num(rd_ests_diff_covs[[x]]$ci[3,2]),
                                   format_num(rd_ests_diff_covs[[x]]$bws[1, 2]),
                                   format_num(rd_ests_diff_covs[[x]]$Nh[1] + rd_ests_diff_covs[[x]]$Nh[2], digits = 0),
                                   format_num(outcome_means[[x]]),
                                   format_num(outcome_std_dev[[x]]),
                                   "differenced covs"))

  names(temp_results) <- c("outcome", "coef", "p_val", "ci_low", "ci_upp",
                           "h", "n", "mean", "sd", "specs")

  results <- rbind(results, temp_results)
}


results_rdrobust <- results

save(results_rdrobust, file = "replication_RD_results_rdrobust.RData")


#####################################
######  FORM OF GOVT ANALYSES  ######
#####################################

## List of DVs
varlist <- c("ttl_revs", "ttl_own_revs", 
             "ttl_taxes", "ttl_sales_taxes", "ttl_prop_taxes",
             "ttl_chgs_miscrev", 
             "ttl_debt", "ttl_debt_issued", "st_debt",
             "ttl_exps", 
             "police", "fire", "admin", "waste", 
             "roads", "parks", "libraries",
             "health", "housing","welfare",
             "basic_services", "developmental", "redistributive")


#################################################
#### RD business exec - mayor-council cities ####
#################################################

mayors_subs <- filter(mayors, !is.na(exec_margin), !is.na(iform), iform == "mayor-council")

results_CCT <- data.frame()

## w/ covs

rdrobust_specs_covs <- lapply(varlist, function(x) {
  
  substitute(with(mayors_subs, rdrobust(i, exec_margin,
                                        covs = cbind(j, 
                                                     log(ipopulation),
                                                     ipctwhite,
                                                     log(imedhhinc),
                                                     log(imedhouseval)))),
             list(i = as.name(paste0("pc_", x, "_lead2")), 
                  j = as.name(paste0("pc_", x, "_lag1"))
             ))
})



rdrobust_results_covs <- lapply(rdrobust_specs_covs, eval)



## these results use CCT optimal bandwidth but regular SEs
for(j in 1:length(varlist)){
  
  ## covs model    
  h_CCT <- as.numeric(rdrobust_results_covs[[j]]$bws[1, 2])
  
  mayors_subs$wgts_CCT <- weights_CCT(mayors_subs$exec_margin, h_CCT)
  
  
  covs <-    lm(substitute(i ~ (exec_margin > 0) * exec_margin +
                             j +
                             log(ipopulation) +
                             ipctwhite +
                             log(imedhhinc) +
                             log(imedhouseval),
                           list(i = as.name(paste0("pc_", varlist[j], "_lead2")), 
                                j = as.name(paste0("pc_", varlist[j], "_lag1"))
                           )), 
                data = mayors_subs, weights = wgts_CCT, subset = !is.na(wgts_CCT) & wgts_CCT > 0)
  
  summ_covs <- summary(covs)
  
  covs_robust <- coeftest(covs, vcov=vcovHC(covs, type="HC1")) 
  
  
  temp_results <- data.frame(cbind(as.character(varlist[j]), 
                                   format_num(summ_covs$coefficients[2,1]),
                                   format_num(covs_robust[2,2]),
                                   format_num(covs_robust[2,4]),
                                   format_num(h_CCT),
                                   format_num(summ_covs$df[2] + summ_covs$df[3], digits = 0),
                                   "CCT covs"))

  names(temp_results) <- c("outcome", "coef","robust_se", "robust_pval", "h", "n", "specs")
  
  
  results_CCT <- rbind(results_CCT, temp_results)
  
}



mayor_council_results_CCT <- 
  results_CCT %>% 
  mutate(form_govt = "mayor_council")



###################################################
#### RD business exec - council-manager cities ####
###################################################

mayors_subs <- filter(mayors, !is.na(exec_margin), !is.na(iform), iform == "council-manager")

results_CCT <- data.frame()


## w/ covs

rdrobust_specs_covs <- lapply(varlist, function(x) {
  
  substitute(with(mayors_subs, rdrobust(i, exec_margin,
                                        covs = cbind(j, 
                                                     log(ipopulation),
                                                     ipctwhite,
                                                     log(imedhhinc),
                                                     log(imedhouseval)))),
             list(i = as.name(paste0("pc_", x, "_lead2")), 
                  j = as.name(paste0("pc_", x, "_lag1"))
             ))
})



rdrobust_results_covs <- lapply(rdrobust_specs_covs, eval)



## these results use CCT optimal bandwidth but regular SEs
for(j in 1:length(varlist)){
  
  ## covs model    
  h_CCT <- as.numeric(rdrobust_results_covs[[j]]$bws[1, 2])
  
  mayors_subs$wgts_CCT <- weights_CCT(mayors_subs$exec_margin, h_CCT)
  
  
  covs <-    lm(substitute(i ~ (exec_margin > 0) * exec_margin +
                             j +
                             log(ipopulation) +
                             ipctwhite +
                             log(imedhhinc) +
                             log(imedhouseval),
                           list(i = as.name(paste0("pc_", varlist[j], "_lead2")), 
                                j = as.name(paste0("pc_", varlist[j], "_lag1"))
                           )), 
                data = mayors_subs, weights = wgts_CCT, subset = !is.na(wgts_CCT) & wgts_CCT > 0)
  
  summ_covs <- summary(covs)
  
  covs_robust <- coeftest(covs, vcov=vcovHC(covs, type="HC1")) 
  

  temp_results <- data.frame(cbind(as.character(varlist[j]), 
                                   format_num(summ_covs$coefficients[2,1]),
                                   format_num(covs_robust[2,2]),
                                   format_num(covs_robust[2,4]),
                                   format_num(h_CCT),
                                   format_num(summ_covs$df[2] + summ_covs$df[3], digits = 0),
                                   "CCT covs"))

  names(temp_results) <- c("outcome", "coef", "robust_se", "robust_pval", "h", "n", "specs")
  
  
  results_CCT <- rbind(results_CCT, temp_results)
  
}



council_manager_results_CCT <- 
  results_CCT %>% 
  mutate(form_govt = "council_manager")


fog_results_CCT <- bind_rows(mayor_council_results_CCT, council_manager_results_CCT)

save(fog_results_CCT, file = "replication_fog_CCT_results.Rdata")








